﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using static TESalgorithm.Constant;
using MathNet.Numerics;
namespace TESalgorithm
{
    public class Equation
    {
        /// <summary>
        /// 普朗克黑体辐射方程
        /// </summary>
        /// <param name="lambda">波长</param>
        /// <param name="T">黑体温度，单位K</param>
        /// <param name="Unit">输入的波长单位,默认是微米</param>
        /// <returns>单位波长辐亮度，单位是W/Sr/m^2/μm</returns>
        public static double Plunk(double lambda, double T, LengthUnit unit= LengthUnit.micons)
        {
            //单位先全部转换成m
            switch (unit)
            {
                case LengthUnit.meter:
                    break;
                case LengthUnit.micons:
                    lambda = lambda * 1e-6;
                    break;
                case LengthUnit.nanometer:
                    lambda = lambda * 1e-9;
                    break;
            }
            //普朗克方程e的系数
            double para1 = h * c / (lambda * k * T);
            double para2 = 2 * Math.PI * h * c * c;
            //普朗克方程计算辐射量度
            double radiance = para2
                              / Math.Pow(lambda, 5) 
                              / (Math.Exp(para1) - 1);
            //转化单位为每微米
            return radiance*1e-6;
        }


        /// <summary>
        /// 普朗克函数的逆函数，根据辐射求温度
        /// </summary>
        /// <param name="lambda">波长</param>
        /// <param name="radiance">辐亮度，单位是W/Sr/m^2/μm</param>
        /// <param name="emmisivity">发射率，如果空缺这一项就默认为黑体</param>
        /// <returns></returns>
        public static double Plunk_inverse(double lambda, double radiance,double emmisivity=1)
        {
            //转化所有单位为每米
            radiance = radiance * 1e6/emmisivity;
            lambda = lambda * 1e-6;
            //逆普朗克函数
            double para1 = h * c / k;
            double para2 = 2 * Math.PI * h * c * c;
            double T_b = para1 / lambda /
                Math.Log(para2 /( radiance * Math.Pow(lambda, 5)) + 1);
            return T_b;
        }

        public static double min_Emmisivity(double MMD)
        {
            return MMD_a - MMD_b * Math.Pow(MMD, MMD_c);
        }

        /// <summary>
        /// 算法的RAT模块，输出五个beita和MMD
        /// </summary>
        /// <param name="emmisivity"></param>
        /// <param name="MMD"></param>
        /// <returns></returns>
        public static List<double> RAT(List<double> emmisivity,out double MMD)
        {
            //平均值
            double average = emmisivity.Average();
            //求每个波段的beita
            var beita = emmisivity.Select(x => x / average).ToList();
            MMD = beita.Max() - beita.Min();
            return beita.ToList();
        }

        /// <summary>
        /// 求出某个亮温、某个波段区间对应的辐射亮度，这里采用了在半波宽度内透过率=1的近似
        /// </summary>
        /// <param name="BBT">亮温</param>
        /// <param name="T">物体的温度</param>
        /// <param name="centre">默认单位为微米</param>
        /// <param name="FWHM">默认单位为微米</param>
        /// <returns>辐射亮度，单位W/m^2/Sr</returns>
        public static double radianceOfBand(double T, double centre, double FWHM)
        {
            double radiance = Integrate.OnClosedInterval(lambda => Plunk(lambda, T),
                centre - FWHM / 2,
                centre + FWHM / 2);
            return radiance;
        }

        /// <summary>
        /// 求出某个亮温、某个波段区间对应的发射率
        /// </summary>
        /// <param name="BBT">亮温</param>
        /// <param name="T">物体的猜测温度</param>
        /// <param name="centre">默认单位为微米</param>
        /// <param name="FWHM">默认单位为微米</param>
        /// <returns>反射率</returns>
        public static double emmisivityFromBBT(double BBT,double T,double centre,double FWHM)
        {
            double BBradiance = radianceOfBand(T, centre, FWHM);
            double radiance = radianceOfBand(BBT, centre, FWHM);
            double emmisity= BBradiance/radiance;
            return Math.Min(emmisity, 0.994);
        }



    }
}
